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Abstract 

We discuss the recently proposed muhicanonical multigrid Monte 
Carlo method and apply it to the scalar (/)^-model on a square lattice. 
To investigate the performance of the new algorithm at the field-driven 
first-order phase transitions between the two ordered phases we care- 
fully analyze the autocorrelations of the Monte Carlo process. Com- 
pared with standard multicanonical simulations a real-time improve- 
ment of about one order of magnitude is established. The interface ten- 
sion between the two ordered phases is extracted from high-statistics 
histograms of the magnetization applying histogram reweighting tech- 
niques. 
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1 Introduction 



First-order phase transitions play an important role in many fields of physics 
[H, 1^ . Examples range from the well-known process of crystal melting through 
the deconfining transition in hot quark-gluon matter to various steps in the 
evolution of the early universe. It is therefore gratifying that recently high 
precision Monte Carlo studies of systems undergoing a first-order phase tran- 
sition have become feasible by showing that the problem of the supercritical 
slowing down, governed by exponentially diverging autocorrelation times 

r oc exp{2aL'^"^}, (1) 

may be eliminated by means of the so-called multicanonical algorithm 
Here a is the interface tension between the coexisting phases, L the linear 
size of a rf-dimensional cubic system, and the factor 2 accounts for the usually 
employed periodic boundary conditions. 

While the multicanonical algorithm does beat the exponential slowing 
down the remaining autocorrelation times typically still diverge with some 
power a ^ 1 ... 1.5 of the lattice volume V = L'^ % |, |, |, |3, |, i, 

r oc V^", (2) 

and may consequently still be severe. It is therefore worthwhile to look for 
further improvements of the Monte Carlo scheme. Critical slowing down with 
a power-law divergence of the autocorrelation time is a notorious problem for 
simulating systems at a second-order phase transition. Fortunately also here 
substantial progress was made in the past few years by designing modified 
Monte Carlo update schemes which reduce or even eliminate the critical slow- 
ing down, see Ref. |T0| for reviews. Among these sophisticated Monte Carlo 



schemes multigrid techniques [11, 12, 13, 14, 15, Qq] have been shown both to 



perform quite successfully and to be rather generally applicable. In a recent 
paper we have therefore proposed |]16|, to combine the multicanonical 



approach with multigrid techniques and have presented preliminary inves- 
tigations ||l^ which show that autocorrelation times in the multicanonical 



simulation may further be reduced by this combination. The purpose of this 
paper is to present a detailed study of the performance of the multicanon- 
ical multigrid method applied to a scalar two-dimensional 0^-theory on the 
lattice. 
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For Potts models it was recently proposed along other lines to combine 
a multicanonical demon algorithm with cluster update methods in a hybrid- 
like fashion |TB[. Another interesting idea is to enhance the performance 



of Monte Carlo simulations of systems at a first-order phase transition by 
treating the parameter which controls the strength of the transition as a 
dynamical variable [|19]. 

The outline of the paper is as follows. After defining the model and dis- 
cussing its basic features in section 2 we will briefly review the characteristic 
properties of multicanonical reweighting and of multigrid update techniques 
in section 3, and show how the two Monte Carlo schemes may be combined. 
We further discuss the error estimates for canonical observables computed 
from multicanonical simulations and introduce an associated effective auto- 
correlation time which allows for a direct comparison with canonical simu- 
lations. Details of the calculation which in principle is straightforward but 
nonetheless requires some care are presented in Appendix A. In section 4 
we focus on the autocorrelation times achieved by the different update algo- 
rithms. After presenting some data for the canonical case for later comparison 
we first analyze autocorrelations in the pure multicanonical distribution. We 
then discuss how the effective autocorrelation time emerges from these data. 
Since the multicanonical multigrid method allows to simulate the model with 
some accuracy, we compute in section 5 the interface tension using histogram 
reweighting techniques. In section 6 we conclude with some remarks on the 
applicability of the method and on future perspectives. 



2 The model and observables 



As a test case to study the performance of the recently proposed multicanon- 
ical multigrid algorithm 
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^ ^ ^ we have taken the two-dimensional scalar 

0^-model on a square L x L-lattice with periodic boundary conditions. This 
model has been studied recently in a number of numerical investigations and 
has repeatedly been used as a testing ground for advanced numerical simu- 
lation techniques ||lT], | 



Previous investigations have focussed 



on properties of the model at criticality ||2^, in particular on the question of 
finite-size scaling and the universality with the Ising model pOl E3[. In this 



paper we present data from simulations performed in the broken symmetry 
phase, that is, strictly below the critical temperature at zero magnetic field 
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(along the first-order transition curve). We focus on tlie autocorrelations 
of the Monte Carlo process and on the interface tension between the two 
ordered phases of positive and negative magnetization. 
The model is defined by the partition function 



exp{-n{M)}, (3) 
where 

nm) = E[^(V0.)^ - ^0,^ + 9<pt]- (4) 
i=i ^ ^ 

Here (V0i)^ = Z])l=i(0i^ — is the squared lattice gradient where de- 
notes the nearest neighbor index in the lattice direction /x, and \^ = L*^ is the 
volume of a o?- dimensional cubic system. The constants fi^,g > are param- 
eters to be specified later, and the energy is measured in temperature units, 
i.e., we have set (3 = l/fc^T = 1. Observables for the model are denoted by 

V 

M, = Y.(Pl k=l,2,..., (5) 

i=l 

and the corresponding densities are denoted by rrii = Mi/V. The kinetic 
term will be denoted by 

Ko ^ E(V0.)^ (6) 

i=l 

Other quantities can be defined as functions of these observables, thus the 
energy E is given by 

E = ^Ko-^M, + gM,, (7) 

and the specific heat C and the (finite lattice) susceptibility x can be obtained 
from 

C = {{E')-{Ef)/V, (8) 
X = {{M',)-{\M,\f)/V. (9) 

In d = 2 dimensions the model (|^,D displays a line of second-order phase 
transitions in the (/i^, 5')-plane which was numerically determined in Ref.[p^ 
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In the thermodynamic hmit the system shows spontaneous symmetry break- 
ing for all /i^ > nl{g) with a nonvanishing expectation value of the aver- 
age magnetization (Mi). Adding a term — to the energy the system 
changes discontinuously from a state of positive magnetization (Mi) > to a 
state of negative magnetization (Mi) < if we tune the magnetic field h from 
positive to negative values. For vanishing magnetic field h = 0, the system 
consequently is at a first-order phase transition. If the linear length is finite, 
L < oo, the system then is flipping back and forth between states of positive 
and negative magnetization which renders (Mi) = also for /i^ > ^1{g). At 
a first-order phase transition point finite systems can also exist in a mixed 
phase configuration which is characterized by regions of the pure phases sep- 
arated by interfaces. For topological reasons there are necessarily an even 
number of interfaces of length L for periodic boundary conditions which we 
have always used. For energetic reasons configurations of more than two in- 
terfaces almost never occur. A typical mixed phase configuration is shown 
in Fig.l. Due to the additional free energy of the interfaces, configurations 
with small total magnetization are suppressed by a factor exp(— 2crL) where 
a is the interface tension. For this reason the probability distribution of 
the magnetization P(mi) shows two distinct peaks separated by a region of 
strongly suppressed mixed phase configurations. 

Following Binder |^ one can extract the interface tension a by determin- 



ing the ratio of the maxima P"^^^ to the minimum p™™ of the distribution 
P{mi). The interface tension is then given by 



with 



a = lim cTi, (10) 



-2a lL 



• (11) 

This expression can easily be evaluated provided that the statistical uncer- 
tainties of both p™^^ and P™™ are small, which can be achieved by perform- 
ing simulations according to the multicanonical distribution. 

Since in this paper we study the 0^-model as a testing ground for the 
performance of Monte Carlo algorithms at first-order phase transitions the 
primary observable of interest will be the average magnetization m\ and 
its autocorrelations in the Monte Carlo process. Although for the partition 
function (H,^ this observable vanishes trivially in finite systems for reasons 
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of symmetry we emphasize that this symmetry is completely artificial and 
would at once be broken, e.g., by adding some term of odd power in to the 
energy in (^. 

Autocorrelations provide a measure for the dynamics of different Monte 
Carlo schemes, see Ref.||2^ for a review. In general, if Oi denotes the i-th 
measurement of an observable O in the Monte Carlo process the autocorre- 
lation function A{j) is defined by 

and from A{j) the integrated autocorrelation time r™* is obtained in the 
large k limit of 

1 k 

r{k) = - + T.Mj)- (13) 
Since for large j the autocorrelation function A{j) decays like an exponential 

A(jV^aexp(-j7r-P), (14) 

where t'^^p denotes the exponential autocorrelation time and a is some con- 
stant, r(/c) behaves like 

oo 

r(k) = r'^'-a ^ exp(-j7r-P) (15) 

j=k+i 

where r™* = r(oo). The latter expression may be used for a numerical 
determination of the exponential and integrated autocorrelation times. Since, 
in general, all these quantities depend on the observable under consideration 
we will indicate the relevant observable by an additional subscript unless it 
is clear from the context which obserable is meant. 

3 Multicanonical Monte Carlo simulations us- 
ing multigrid techniques 
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3.1 Multicanonical simulations 



The basic idea of the multicanonical approach p, ^ |^ is to simulate an 
auxiliary distribution in which the mixed phase configurations have the same 
weight as the pure phases and canonical expectation values are recovered 
by reweighting. Hence the multicanonical approach is not itself a Monte 
Carlo update algorithm but a general reweighting prescription which allows 
to simulate distributions which are numerically easier to handle. 

While similar ideas have been known in the literature under the name 



of umbrella sampling already for a long time p6[ the practical relevance 
of multicanonical reweighting techniques for simulations of first-order phase 
transitions ||, |, ||, |^ was realized only recently 0. The multicanon- 
ical approach may be formulated in a rather general way for a variety of 
applications but in the context of first-order phase transitions and quan- 
tum mechanical tunneling problems [p!6| , |1^] it may simply be regarded as 
being basically a reweighting technique 0. 

Let m = m{{(j)i}) be an observable whose probability distribution in the 
canonical ensemble displays two strong separated peaks. In a field driven 
transition, as in our case, m({0j}) is the magnetization mi and the situation 
has also been referred to as multimagnetical simulation ^j. For a temperature 
driven transition the relevant observable would be the energy of the system. 
In the multicanonical reweighting approach we now rewrite the partition 
function by introducing some function / as 

Z = Y[ [ t/0.e-««'^'»-/(-)e^('"), (17) 
i=i 

and adjust the reweighting factor exp(— /(m)) in such a way that the re- 
sulting histogram of m sampled according to the multicanonical probability 
distribution 

p~{M) oc exp[-7^({0,,})-/(m)] 

is approximately fiat. Here Ti™"*^** is the central object of a multicanonical 
simulation, and plays the same role in it as Ti. does in a canonical simulation. 
Canonical observables (O)can can be recovered according to 

(O)... ^ ^. (19) 
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where (...) without subscript denote expectation values in the multicanonical 
distribution and w = exp(/(m)) is the inverse reweighting factor. 

The multicanonical probability distribution p™^'^^ may be updated us- 
ing any legitimate Monte Carlo algorithm. The simplest choice is a local 
Metropolis update where we consider as usual local moves (pig + A0jg 

at some site io and compute the energy difference Ai?™"'^^ according to (p!8|), 
i.e., the decision on whether a Metropolis move will be accepted or not is 
now to be based on the energy difference 

AE~ = AE + f{m + Am) - f{m), (20) 

where AE = ?^(0i, 02, • • • , <Pio + ^0*0; • • • ; ^v) — '^({0t}) is the canonical 
energy difference and Am = m(0i, 02, • • • , 0io + A0io, • • • , 'Pv) — '^({0j}) is 
the corresponding difference in the observable m. If the canonical probability 
distribution is reweighted in the magnetization, m({0j}) = mi = Y.Y=i 0j/^, 
this difference is simply given by Am = Aipi^/V . For a temperature-driven 
transition we have m = Ti/V and Am = AE/V. 

In practical simulations / may be recorded in form of a simple stair case 
function which does not introduce any numerical inaccuracy since this factor 
cancels out in all canonical expectation values. It is also worth mentioning 
that since m depends on all values of 0i, the resulting multicanonical energy 
is essentially nonlocal. 

As we will discuss in Section 3.3, the multicanonical probability distribu- 
tion p'^™^''^ may also be updated by a multigrid Monte Carlo method. 



3.2 Multigrid techniques 

The basic idea of multigrid Monte Carlo techniques [0, |T2|, [1^, |16 



2% P8| is to systematically perform updates on different length scales of the 
system. In the corresponding unighd viewpoint, which always looks at the 
effects on the original fine grained lattice, this is done by moving blocks of 
1, 2"^, 4'^, 8'^, . . ., 2"'^ = V adjacent variables at a time. In the multigrid 
formulation these collective update moves are implemented by introducing 
auxiliary fields on coarse-grained lattices. Specifically one introduces a se- 
quence of coarsened grids = n — 1, . . . , of size 2'^''. In the simplest 
piecewise constant interpolation scheme we identify a pair, square, cube, etc. 
of 2'^ neighbouring grid points on a grid of level k with a single grid point of 
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the next coarsened grid On these coarsened grids we have auxihary 

fields 0^'^"^'* representing the collective moves of the original fine-grained lat- 
tice. A (piecewise constant) interpolation operator V is defined by simply 
adding a finite value of some variable (pf^'^^ on a coarse grid to each of the 
2d corresponding grid points of the next finer grid. This allows to define a 
Hamiltonian of the coarse grid in terms of the Hamiltonian on the next finer 
grid by [|lT 



n^'-'\{<Pt'^}) = ^^'^{#1 + n{<pt'^})). (21) 

In essence this prescription defines a Hamiltonian on the coarse grid H^'^"^^ 
by freezing the field variables of the next finer grid and calculating the 
effect of collective moves represented by the field variables added onto 

S'^'^^ by the piecewise constant interpolation operator. If the functional form 
of the Hamiltonian remains stable under the coarsening prescription 
this multignd implementation of the collective move updates minimizes the 
amount of computational effort compared to the straightforward unigrid im- 
plementation of the collective move update, in a way similar to the Fast 
Fourier Transformation (FFT). Also, it allows to define the multigrid update 
in the following recursive way. Updates of level S^'^^ consist of a) rii presweeps 
using any valid local update scheme with Hamiltonian (|2l|), b) calculating 
the Hamiltonian for the next coarser grid H*^'^"^^ (which according to (pi|) 
depends on the current configuration on grid S*-'^'-') and initializing the vari- 
ables on grid S^^^^^ to zero. One then c) updates the field variables 4>[''~^^ 
by applying the multigrid update ■jk-i times. To complete the update cycle 
one then d) interpolates the variables of grid S^^~^)back to grid S*^*^-' and e) 
performs another n2 postsweeps of the local update algorithm. On the coars- 
est grid, of course, one only performs steps a) and e). In this way we cycle 
through the sequence of coarsened grids in a specific manner which is deter- 
mined by the parameters jk- Particularly successful is the choice 7^ = 2, a 
sequence which is commonly called W-cycle since its graphical representation 
very much resembles the letter W.Q 



iSee, e.g., Ref.M, p.33. 
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3.3 Multicanonical multigrid Monte Carlo 

From the unigrid viewpoint it is immediately clear that the multicanonical 
and multigrid methods can easily be combined for a field-driven transition 
where m = mi is the average magnetization [T^. Since on level k we 
effectively always move 2^"'~''^^ spins in conjunction an accepted Metropolis 
move would change the average field mi by an amount of 2^''-'''^'^A(j)[''J/V. 
The only modification for the update on level k will therefore be to compute 
the energy difference according to 

^^muca,{fc) _ ^^{k) + + 2("-'=)'^A0jVv^) - /(m), (22) 

where AE~^^''^ is the energy difference computed with the coarse-grid 
Hamiltonian (^) as in the usual canonical multigrid formulation. While this 
modification is obvious from the unigrid point of view it should be stressed 
that the modifications for a recursive multigrid implementation are precisely 
the same. 

In our case the doubly peaked observable m which controls the multi- 
canonical reweighting factor is the magnetization mi. In this case the nec- 
essary modifications for a recursive multigrid update of the multicanonical 
Hamiltonian Ti™"'^'^ are in fact almost trivial. The combination of multigrid 
update schemes with the multicanonical reweighting idea, however, is nei- 
ther restricted to the special choice of m = mi nor to any special choice of 
the reweighting factor /. In general, a multigrid Monte Carlo update of a 
multicanonical Hamiltonian H^^^^{{(j)i}) = + should be 

feasible and effective whenever both the canonical Hamiltonian Ti and the 
parameter m are stable under coarsening. To see this, let us assume we would 
want to reweight the canonical Hamiltonian Ti. in some other observable of 
{(pi}, say m = m2, rather than in m = mi. In order to compute the reweight- 
ing factor f{m2) on level k we would need to know the actual value of m2 as a 
function of the coarse-grid variables </)|^\ Clearly, we cannot simply compute 
the average value of m2 by multiplying A^j^^^ with a simple factor as was the 
case for the average magnetization. Indeed, from the unigrid point of view 
we would need to know the actual configuration of the original-grid variables 
and the efficiency gained from the multigrid update would be lost at least for 
this part of the update. It is therefore important to realize that m2 may also 
be calculated using the usual coarsening prescription. In general, the analog 
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of eq. (^) for the function m would simply read 

^(.-1) (1^(^-1) I) ^ ^W({0f )} + (23) 

Therefore we would only have to compute the coarse-grid coefficients for the 
function m'-'^"^^ in addition to those for and we could then use the 

value of m^^'^^ in order to compute the reweighting factor. From this point of 
view, the factors 2^'^^'^^^ appearing in eq. ( P^] ) are nothing but the coarse-grid 
coefficients for the average magnetization using piecewise constant interpo- 
lation. 

We would like to stress again at this point that the condition that m 
remains stable under coarsening is the only restriction for an effective multi- 
grid Monte Carlo update of a multicanonical Hamiltonian. In particular, this 
means that a) the reweighting factor f{m) may be computed for any func- 
tion /. In fact, the step functions normally employed are examples for rather 
special, highly non-linear functions, b) The condition also allows a combi- 
nation of multigrid techniques and multicanonical reweighting if we take the 
canonical Hamiltonian Ti itself as the observable m. This would be the sit- 
uation for a temperature-driven first-order transition. If Ti. is stable under 
coarsening (as we always assume it is) it therefore immediately follows that a 
multicanonical multigrid Monte Carlo simulation would be perfectly feasible 
in this case as well, c) We also believe, that multicanonical multigrid Monte 
Carlo simulations should be feasible for models other than those characterized 
by a Hamiltonian of the form (^) such as non-linear 0(n) sigma- models. If, 
e.g., one would want, for some reason, to simulate an XF-model with a mul- 
ticanonical reweighting factor / = f{sx) where Sx is the average x-component 
of the spins one might express Sx as cos{Q ij)/V and the latter function 
is stable under coarsening for piecewise constant interpolation if we allow 
for an additional term (I]iSin(0j))/\^ on the coarsened grids. More realistic 
but nevertheless feasible as well would be the case where the reweighting 
variable is the squared magnetization d) Furthermore, the fact that we 
only need to know the actual value of m for each coarse-grid update and 
the fact that this value may be computed from eq. ( PBj ) also entails that the 
multicanonical multigrid Monte Carlo method is, in general, not restricted 
to the piecewise constant interpolation scheme, e) Finally, it should be re- 
called that multigrid techniques are sophisticated update schemes which do 
not presuppose specific update algorithms. In principle, we may therefore 
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use any other valid update algorithm on the coarsened grids instead of the 
Metropolis update. 

In this paper we will substantiate our claim that multicanonical multigrid 
Monte Carlo simulations are both feasible and profitable by a careful analysis 
of the performance of the method for the model @, (Q). For other situations 
the method should be tested by explicit simulations. As long as the central 
condition for the multicanonical multigrid approach is fulfilled, however, we 
do not expect any difficulties regarding the feasibility of the method. 



3.4 Effective autocorrelation time 

Before discussing our results it is worthwhile to comment on a technical 
complication which arises in evaluating the efficiency of multicanonical sim- 
ulations. In previous investigations it was the exponential autocorrelation 
time measured in the multicanonical distribution which was used to estimate 
the performance of the multicanonical algorithm. Alternatively, autocorrela- 
tions were also measured by counting the average number of sweeps needed 
to travel from one peak maximum to the other and back (see section 4.1. 
below). This nicely illustrated the absence of an exponential slowing down 
[0, 1^, 1^, ^ D, ^ . It should be stressed, however, that neither the exponen- 
tial autocorrelation time nor the diffusion timeQ can a priori serve as a fair 
quantitative measure for comparison of the performance of multicanonical 
simulations with canonical update schemes. The reason is that the estimator 
O = X^il™ OiWi/ Yl!i=i_Wi of canonical observables (|19D is a ratio of two dif- 
ferent multicanonical observables which may have different autocorrelations 
and, moreover, are usually strongly cross-correlated. It is thus not immedi- 
ately obvious how autocorrelations relevant for canonical quantities should be 
defined and measured in multicanonical simulations. For a fair comparison 



with canonical simulations we therefore define an effective autocorrela- 
tion time t"^^ and write the error estimate also in the multicanonical case in 
the standard form 

^6 = (^oJ (24) 



^In analogy to canonical simulations this is often called "tunneling time" even though 
this is quite a misleading terminology in the multicanonical case where the dynamics is of 
a diffusive type. 
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where Nm is the number of multicanonical measurements and (a^J^"" is the 



variance of Oi in the canonical ensemble computed according to eq. (|T9|) . A 
simple and numerically stable way to obtain an estimate for the effective 
autocorrelation time is to compute the error of the estimator O by 
jackknife blocking. In principle, r^f can also be calculated by applying stan- 
dard error propagation starting from the basic reweighting formula (plQl). As 
shown in Appendix A the squared canonical error of the canonical esti- 
mator O of an observable (O)can based on multicanonical measurements 
is then given by 

\'-^/can [ //n , . \9 at 



(25) 



where {x; y) = (xy) — (x) (y) is the covariance matrix of two observables x and 
y. From this expression it is clear that in general three different integrated 
autocorrelation times Tq^.q^, t™1„, and Tq^.^ of the observables Ow and w 
(measured in the multicanonical distribution) must be taken into account. 



4 Results 

We have simulated the model (^^ in two dimensions {d = 2) at three differ- 
ent points in the (/i^, 5f)-plane, namely we have taken g = 0.25 and varied /x^ 
as /i^ = 1.30, 1.35, and 1.40. For this value of g, the second-order phase tran- 
sition to the disordered phase occurs at /x^ = 1.265(5) as it was determined 
in Ref. ||2^, confirmed in Ref. [^, and reproduced in our own simulation 



(see section 5). With this choice of parameters our simulations were per- 
formed in a regime which shows the typical first-order behavior already for 
relatively small lattices. For large linear lattice size L the ratio between the 
maxima and the minima of the histogram for mi will then easily take on 
several orders of magnitude. For the severest case which we have investi- 
gated (/i^ = 1.40, L = 64) this ratio, e.g., is already more than 9 orders of 
magnitude. In these extreme cases an important point of the multicanonical 
algorithm is the way to obtain the trial histogram since the performance of 
the multicanonical simulation strongly depends on the quality of the assumed 
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trial distribution. If there is no chance to make an initial guess on the basis 
of some knowledge of the system, the most straightforward way to proceed is 
by iterations which in our case was done as follows. Starting with a canonical 
simulation we performed some thermalizing sweeps and then obtained a first 
histogram on the basis of 50 000 x ng sweeps. Here Ue is a parameter which 
allows to adjust the time scale of the MC process, i.e., measurements were 
always taken only after every "empty" sweeps. In the multigrid case we 
count a complete cycle as one sweep, and we only performed presweeps, i.e., 
we always had ni = 1 and n2 = 0. In order to maintain a roughly constant 
Metropolis acceptance rate of about 50% we had to scale down the maximal 
step width oc 0.6""^^ which conforms with recent analytical investigations 
of the Metropolis acceptance rate [PT|. The first histogram was then sym- 



metrized and any empty bins between the peaks were filled by interpolation 
using rough estimates for the interface tension obtained from simulations 
of the smaller lattices. Using this histogram as a first guess to construct 
the reweighting factor exp(— /) we performed another 50 000 x Ue multi- 
canonical sweeps. The resulting histogram proved to be sufficient for lattice 
sizes L = 8, 16, and 32. To obtain high precision the resulting histogram 
nevertheless was in any case again symmetrized and taken for the final sim- 
ulation run. For L = 64 we have done one more iteration of 1 000 000 x Ue 
sweeps and used only this resulting histogram as trial distribution for the 
final run. For the determination of the trial histograms, which always had 
a bin size of 0.008, we have in any case used the W-cycle and, to allow for 
a direct comparison, we have taken the same trial histograms for the final 
runs of both the standard multicanonical simulation and the combination 
with the multigrid scheme. In our final runs we have in each case performed 
1 000 000 X He sweeps after discarding 10 000 x ne sweeps for thermalization. 
To allow for later histogramming and fiexible analysis of autocorrelations we 
have recorded the time series of Kq, Mi, M2, and M4, and all errors were 
computed by jackkniving the data with 50 blocks. 

Figures 2(a-c) show the fiat multicanonical distributions and the canoni- 
cal double-peak histograms of mi after reweighting for /z^ = 1.30 and different 
lattice sizes L = 8, 16, and 32. The quality of the flatness of the multicanoni- 
cal distributions for our largest lattices (L = 64) and for /i^ = 1.30, 1.35, and 
1.40 can be judged from Figs. 3(a-c). Although the multicanonical distribu- 
tions are essentially flat between the peaks there still is some structure in the 
multicanonical distributions which affects autocorrelations of mi. Note the 
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flat regions in the canonical histograms for mi ^ which reflect the fact that 
mixed phase configuration for some range of mi all look similar to the one 
displayed in Fig. 1. Different total magnetizations in this regime result from 
a relative shift of the two interfaces and the free energy of these configura- 
tions is the same as long as interactions between the interfaces are negligible. 
The arrow in Fig. 3c indicates the value of mi for the configuration displayed 
in Fig. 1. 

The drastic difference between the canonical and the multicanonical up- 
dates can be illustrated by looking at the time series of mi as shown in 
Figs. 4 (a-d). While the time series in the canonical case clearly displays 
the instantons characteristic for tunnelling processes the observable in the 
multicanonical simulation shows a random walk like behavior. In either case 
the use of the multigrid update affects the time scale of the autocorrelations. 
Note that the time scale in the figures has been adjusted in such a way that 
they display the time evolution over a length of roughly 30 x r™* in either 
case (cp. Tables |] and |^ below). 

To get a more precise view of the performance of the different Monte 
Carlo schemes we have measured autocorrelation times of the relevant ob- 
servables. In order to obtain estimates for the integrated autocorrelation 
time r"* a common way to proceed is to cut the sum ( |13|) self-consistently 
at k = ricnt X where ncut is usually chosen to be 6 or 8. As long as 
this method usually gives sufficiently reliable values. If, however, 
r^'^P is appreciably larger than r™* this method systematically underesti- 
mates the integrated autocorrelation time. Let us illustrate the problem for 
the rather extreme case of (9 = mi exp(/), /i^ = 1.40, and L = 64, cp. Fig.5. 
Here we find an exponential autocorrelation time r^^^ = 3330(530) (in units 
of cycles) which is more than 4 times larger than the true integrated autocor- 
relation time r'°* = 778(63) (cp. Table |^ below). Computing the integrated 
autocorrelation time by self-consistently cutting the sum in eq. ([T^) therefore 
underestimates r™* to be 505(15) for ricut = 6 and 627(25) for ricut = 8. In 
order to circumvent this problem we have therefore proceeded as follows (see 
Fig. 5). For the exponential autocorrelation time r^^^ we have first obtained 
a rough guess t'^^P'^^^ by a linear fit of \nA{j) from j = r™* . . . 3r"^* where 
r"^* is the integrated autocorrelation time obtained by cutting eq.([T3|) self- 
consistently with rZcut = 8- The inset in Fig. 5 shows a logarithmic plot of the 
autocorrelation function A{j) for the example discussed above together with 
this first rough guess shown by the dashed line. Clearly A{j) does not behave 
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like a single exponential as would be the case if r™* = t^^^ and consequently 
one must be careful to fit In A{j) sufficiently far away from the origin. For 
our first approximation we obtained in this case t'^^P'^^-' = 2480. We have 
then performed a three-parameter fit of T{k) according to eq.(|16D in the range 
k = r^^P'W . . . ^t-<^^pXo) -which yielded the values for r'°* and t'^^p quoted in 
Tables 1^-^. Fig. 5 shows T{k) and the corresponding fit. The horizontal lines 
represent our value of r™* together with its error bounds and one clearly sees 
that r(/c) saturates towards this values for k oo. Again we see that the 
data do not yet saturate for k = 6200 ~ 8r'°*. The solid line in the inset 
shows a linear fit of In A{j) in the same range j = t^^P'^^^ . . . ^t^^pM^ This fit 
yields a value r*^^P = 2780(670) which is consistent with the one we quote but 
has larger error bounds. Summing up, we note that fitting r(/c) according to 
eq.([T6|) rather than fitting In A{j) produces simultaneously unbiased values 
for both t'^^p and r™* with smaller statistical uncertainties. Also these fits 
were satisfactorily stable against variations of the fitting range. As usual, 
error bars for the values of r™* and t^^^ obtained by this fitting procedure 
were obtained using the jackknife method. 

4.1 Autocorrelations in canonical simulation 

For later comparison with the multicanonical case we have first performed 
a number of canonical simulations using both the standard local single-hit 
Metropolis update and the multigrid W-cycle. Since the main focus of this 
paper, however, was to investigate the improvement gained by combining 
the multicanonical approach with multigrid techniques, standard canonical 
simulations were done only for /i^ = 1.30 and only for small lattices [L = 
4,8,16). Table ^ shows the measured autocorrelation times for O = mi 
and 777-2 (always given in units of sweeps resp. cycles). For O = mi we see 
that within the error bars the integrated autocorrelation times do not differ 
from the exponential autocorrelation times, i.e., in this case we are dealing 
with an almost purely exponential autocorrelation function (see below for 
a theoretical explanation of this behavior). For the even observable O = 
7772 we find on the other hand that the exponential autocorrelation time 
t'^^p is appreciably larger than the integrated autocorrelation time r™*. The 
difference between r™* and t'^^p increases on larger lattices. Comparing the 
improvement of using multigrid techniques with a W-cycle gives a factor of 
roughly 10 ~ 20. An improvement of this order had already been found in 
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Ref. |Tl|] for simulations at the critical line. In either case the autocorrelation 
times, however, diverge exponentially with increasing linear lattice size L. 

To obtain a rough estimate of the order of magnitude of the quantities 
involved we have fitted the integrated autocorrelation times of nii according 
to r'^* = const x L°exp(2aL). We find t'""^ = 6.41 x L^-^^ exp{2 x 0.027L) 
for the Metropohs case and r''^* = 4.83 x L°-24exp(2 x O.IOL) for the W- 
cycle. Since we expect that the exponent should depend on the interface 
tension we have also performed a constrained fit using for a the value for the 
interface tension obtained from our multicanonical simulations (cr = 0.03443, 
see below). For this fit we find t'""^ = 7.30(36) x L^-^^^^^') exp{2aL) for the 
Metropolis case and r''^* = 1.638(52) x L^-^^eCie) exp(2aL) for the W-cycle. 
Clearly, these fits can only give a rough estimate of the magnitude of the 
relevant quantities for two reasons. First we have only three data points 
for fitting two resp. three parameters. Second, we expect corrections to be 
still appreciable at least for the smallest lattice used (L = 4). Nevertheless, 
comparable numbers were found for the two-dimensional Ising model where in 
[^] a behavior of r™* = 6.80L^'^^ exp(2 x 0.185 x L) was found for the canonical 
heatbath. For the 7-state Potts model fits of the form I.OIL^'^^ exp(2 x 
0.01174 X L) have been reported in Ref. 



The dynamical origin of the autocorrelation time for mi may be illus- 
trated by a simple two-state flip model. We measured the mean time of 
staying in one of the potential wells by digitalizing the evolution series cor- 
responding to a simple two-state model with single flip dynamics This 



procedure is illustrated for the time series in Figs. 4a and 4c. Counting the 
number of Monte Carlo time steps that the systems needs to flip from one 
state to the other we measure an average flip rate. A theoretical analysis 
of this single flip dynamics shows that the exponential autocorrelation time 
for this model is given by t^'p where Ar^^"^ is the average time the system 
needs to flip from one state to the other and back. The measured flipping 
times are also shown in Table For large t^^p the variance of t^^p is 
given by (r*^^P)^ itself as can be calculated in the single flip dynamics and as 
we have verified in our canonical simulations. The error bars for the flipping 
times r*^'P therefore were calculated as 6t^^^ = r^^^/y/n^ where l/nfjip is the 
measured total flip rate, i.e. 2T^^^nf{ip = Nm x Ue- As can be seen in Table 
m application of the multigrid algorithm speeds up the Monte Carlo process 
by increasing the average flip rate by a roughly constant factor. 
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4.2 Autocorrelations in multicanonical simulation 

In the multicanonical simulation the exponential supercritical slowing down 
is overcome by simulating an auxiliary distribution which is flat between 
the two maxima of the histogram, cp. Figs. 2 and 3, and in this case we 
expect the autocorrelations to be governed by some random walk dynam- 
ics, cp. Figs. 4b and 4d. Since multigrid techniques can be applied in the 
multicanonical distribution as well it is of interest to see whether a further 
reduction of autocorrelations can be achieved by this combination. From 
the rather different scales of Figs. 4b and 4d it is already clear qualitatively 
that autocorrelations are in fact reduced. In order to see quantitatively how 
multicanonical simulations are improved by multigrid updating we first mea- 
sured autocorrelations of the corresponding observables in the multicanonical 
distribution using both the standard single-hit Metropolis algorithm and the 
W-cycle. Table |^ shows our results for this case. We see that for O = mi 
it is again found that r™|^ ^ '^mi^ i-^-' ^^so in the multicanonical dynamics 
the autocorrelation function decays like a pure exponential. For the even 
observable O = m2, on the other hand, there is again a difference between 
integrated and exponential autocorrelation times. 

Comparing the absolute values for the Metropolis update and for the W- 
cycle we find that for both observables the multigrid method reduces the 
autocorrelation times by a factor of roughly 15 ~ 20. 

We have also looked at the lattice size dependence of the autocorrela- 
tion times by fitting the data for O = mi according to r™* = aint-Z^^"'* or 
r°^P = Oexp-Z^^'^'P (where z = ad). Here we first note that trivially the expo- 
nent Zint for the integrated autocorrelation times agrees with the exponent 
Zexp for the exponential autocorrelation times. We also find that fitting only 
the data for L = 8, 16, and 32 yields approximately the same exponents for 
the Metropolis case and for the multigrid W-cycle. Looking at the depen- 
dence on the parameter /x^ we find that the exponent increases with /i^, i.e. 
with the strength of the transition. Fitting the integrated autocorrelation 
times obtained from the multigrid simulation for lattice sizes, L = 16, 32, 
and 64, we obtain exponents of about 2.2, 2.5, and 2.7 for /i^ = 1.30, 1.35, 
and 1.40. Including the L = 8-data worsens the fits, and we therefore believe 
that for a reliable estimate one would need to include even larger lattices. In 
general, we observe that the fits have rather large chi-squares, and we hesi- 
tate to draw any definite conclusions. The deviations from linearity found in 
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the log-log-fits are believed to be effected by the fact that the multicanonical 
histograms are not ideally flat. In Figs. 2 (a-c) and 3 (a-c) the multicanonical 
histograms look better than they actually result of the logarithmic 

scale. On a linear scale one still discerns some structure in the supposedly 
flat region between the peaks at least for the larger lattices. Therefore the 
statistical accuracy of our data seems to be better than the systematic fluctu- 
ations of the autocorrelation times caused by imperfect multicanonical trial 
histograms. 

With respect to the purely multicanonical dynamics, we conclude that 
the multigrid technique does not affect the exponent z = ad. However, 
it is by largely reducing the overall scale of the autocorrelation, i.e. by 
reducing the prefactor a, that application of multigrid techniques gives an 
improvement factor of roughly 15 ~ 20, i.e. of about one order of magnitude. 
The improvement factor shows a weak tendency to increase if /i^ approaches 
the critical value. 



4.3 Effective canonical autocorrelations in multicanon- 
ical simulation 

Clearly, the multicanonical reweighting factor is an algorithmical artefact 
introduced in order to obtain higher statistical accuracy for the measurement 
of canonical observables. For a fair comparison between the canonical and 
the multicanonical simulation we therefore have to estimate the error bars 
associated with the canonical observables. 



Odd observables: For odd observables standard error propagation start- 
ing from eq.(19) shows that the effective autocorrelation time for canon- 



ical estimates obtained by measurements in the multicanonical distribution 
is given by 

= J^^j^-TOw^Ow, (26) 

(see Appendix A) with the effective multicanonical variance 
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Table |^ shows the measured integrated and exponential autocorrelation times 
r"^* and t'^^p for O = mi exp(/). Also shown are the effective multicanonical 
variance cr^ according to eg. (P7|) , the canonical variance {o'oj'^'^^ of O = mi, 
which can be computed in a multicanonical simulation by using eg. (p!9D , the 
effective autocorrelation time computed according to eg.(^), and the 
diffusion time r^f defined in analogy to section 4.1. 

First, we notice that, within error bounds, the exponential autocorrelation 
times for O = miexp(/) agree with the purely multicanonical autocorrela- 
tion times for O = mi listed in Table |[ This is not surprising since we 
are still dealing with an odd observable whose slowest autocorrelation mode 
should be same. The integrated autocorrelation times on the other hand 
in this case differ appreciably. This is an indication that the autocorrela- 
tion function (|l^) does not behave like a simple exponential oc exp(— j/r'^^P). 
Rather it is composed of many different modes with only the slowest mode de- 
caying with r^^P as illustrated in the inset of Fig. 5. The relative difference 
between the integrated and the exponential autocorrelation time increases 
both with the size of the system and with /x^. The ratio does not depend on 
the other hand on the use of the algorithm, being roughly the same both for 
the standard Metropolis update and for the multigrid update. 

Table ^ lists both the effective multicanonical and the canonical variances. 
These allow to compute the final effective autocorrelation times which are 
also reported. While the canonical variance depends only weakly on the size 
of the system the effective multicanonical variance cx^ varies appreciably with 
the linear lattice size L. In the worst case, fi^ = 1.40 and L = 64 the ratio 
is already cr"^/ {ctq.Y^^ ~ 11. Conseguently the effective autocorrelation time 
which should be used for comparisons with canonical algorithms is much 
larger than the simple exponential autocorrelation time t'^^p. 

To allow for further comparison with the literature, we have looked also 
for this situation at the exponents Zint resp. Ze^p of the power-like divergence 
r = aL^. In contrast to the purely multicanonical case the exponents Z[^t here 
differ from the exponents Zg^p. While the exponents Zexp for O = miexp(/) 
roughly agree with those for the purely multicanonical observable O = mi 
and thus increase with /x^, the exponents Zi^t seem to stay constant with in- 
creasing /i^. Finally we note that the exponent z^s for the effective autocor- 
relation time r*^^ again increases with //^ which directly reflects the scaling of 
the ratio a^/ {cto Y^^ with /i^ and L. The increase of z^ii with /i^ as compared 
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to Zint is illustrated in Fig. 6. Here the integrated autocorrelation times r"* 
for O = mi exp(/) are shown together with the corresponding effective auto- 
correlation times Teff. The straight lines show fits of the form r™* = Oint-Z^^'"* 
and T^^ = aeffL^<=ff. From this figure it can also be seen that the data for 
the smallest lattice, L = 8, need still be excluded to obtain satisfactory fits. 
Fitting the data for L = 16, 32, and 64 we find for the effective exponents Zes 
values of about 2.3, 2.7, and 3.0 for /x^ = 1.30, 1.35, and 1.40. The exponents 
obtained in our work confirm qualitatively the exponents found for standard 
multicanonical simulations of the two-dimensional g-state Potts model where 
exponents of z = 2.65(5) for g = 7 [^] and z = 2.65(2) for g = 10 have 
been obtained from analyses of the diffusion times. Note that for a random 
walk like behavior as in the multicanonical case one cannot unambiguously 
identify distinct states anymore. One often employed possibility is to mea- 
sure the average number of multicanonical sweeps or multigrid cycles needed 
to travel from one (canonical) peak maximum to the other and back. In 
analogy to the definition for canonical simulations the r^f given in Table ^ 
are one quarter of this average travel time. By using this definition of we 
obtain a nice agreement with at least for the large lattices. A priori, how- 
ever, other definitions of t'^'p are reasonable as well (e.g., using (|mi|)can for 
the cuts instead of the peak locations), and it is difficult to argue which one 
should give the best quantitative agreement with the unambiguously defined 
effective autocorrelation time t'^^. For this reason a direct measurement of 
t'^^ is to be preferred rather than any analogue of the two-state flip model. 

For odd observables the distinction between the directly measured inte- 
grated autocorrelation time and the effective autocorrelation time does not 
pertain to the comparison between the standard multicanonical Metropolis 
update and combination of the multicanonical approach with multigrid tech- 
niques. Since the multicanonical approach is a mere reweighting technique 
0"^ and {o'OiY'^^ affected by applying different update algorithms. 

Hence we find indeed that the same improvement factors of about 15 ~ 20 
are gained both for the integrated autocorrelations and for the effective auto- 
correlation times.0 These effective improvement factors are slightly smaller 
than those found for O = mi from Table and again show a tendency to 
increase when /i^ approaches the critical value. 

Apart form work estimates to be discussed below. 
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Even observables: For even observables we cannot exploit the fact 
that (OiWi) vanishes identically for reasons of symmetry in order to simplify 
the error propagation formula (pS]). In general, to obtain estimates for the 
canonical error of an observable O we therefore have to take recourse to the 
full expression of error propagation (^). Table | shows our results for the 
even observable O = m2 from our simulations at /i^ = 1.30. Here we list 
measured values for all quantities which enter the error propagation formula 
(pSl). The mean values, variances, and covariance, for the Metropolis case 
and for the W-cycle, are consistent within error bars, as, of course, should be 
the case since these quantities do not depend on the update algorithm. The 
integrated autocorrelation times, however, again differ by a factor of roughly 
20. We did not list the corresponding exponential autocorrelation times since 
these agree, within error bounds, with the exponential autocorrelation times 
for O — 7712 listed in Table ^. We have also checked that Tq^.^ ~ '''urOw 
would be expected because of time reversal invariance. 

Next to these values we then list in Table § the squared statistical error 
for the canonical estimator of O = m2 calculated by error propagation from 
the data listed before. Note that the error of this error, however, as well 
as all the errors given in the Table were not computed by error propagation 
but, as usual, calculated directly by jackkniving. 

Clearly, it is this statistical error for the canonical estimates of the ob- 
servable which one wants to reduce by sophisticated Monte Carlo methods. 
When interpreting the statistical errors reported in Table § the time scale 
set by He should also be taken into consideration. While the squared errors 

are approximately of the same order for the Metropolis (M) update and 
for the multigrid W-cycle update (W) we also had to perform many more 
Metropolis updates since we had adjusted ng. If, e.g., for L = 8 the statisti- 
cal error for the Metropolis update is only twice as large as the one for the 
multigrid update, we also had rie = 5(M) resp. 1(W), cp. Table |[ Therefore 
the improvement is given by (0.540/0.2442) x 5 ~ 11 which is roughly of the 
same size as the ratio of the measured autocorrelation times. 

Another technical remark is due at this point. Applying formula (^) 
to calculate the canonical error of multicanonical measurements we run into 
a nasty problem of numerical cancellation. To illustrate this cancellation 
problem let us look at the data for L = 64. Here we find for the first two 
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terms in eg. (p5|) 
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and for the third term 
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i.e., we have a numerical cancellation up to the fifth digit. Also if we look 
at 0"^ we find the same problem which should not come as a surprise if we 
recall that in the definition ( |A10| ) of o"^ we simply dropped the r's in eq. (|25|) . 
Since from Table | we see that for the even observable m2 we always have 
'Tqw-Ow ~ '^w^w ~ '^Ow-w the numerical cancellation should therefore carry over 
to 0"^ as well. Consequently, we obtain numerical results for the statistical 
error estimate which may be completely erroneous. In fact, for L = 64 the 
effective autocorrelation time turns out to be negative which, of course, is 
complete boloney. Therefore it is somewhat difficult to judge the quality 
of the performance of the multicanonical simulation of even canonical ob- 
servables by applying error propagation. Alternatively we can, of course, 
judge the improvement gained by applying multigrid techniques by compar- 
ing the errors obtained by jack-knife blocking procedures. For comparison 
we therefore have listed the squared canonical errors obtained in this manner 
as well as the effective autocorrelation times derived from these jackknife 
errors. These values in general turn out to agree roughly with the calculated 
errors for small lattices but deviate strongly for our large lattices. In general 
we tend to believe that in this case the error estimates obtained by direct 
jackkniving are more reliable than the ones calculated by error propagation. 

Finally it should be remarked that a measurement of even observables in 
the multicanonical distribution is somewhat academic anyway since they may 
already be measured quite accurately in the canonical distribution. Compar- 
ing the autocorrelation times given in Table |^ with the autocorrelation times 
for the canonical simulation reported in Table |I] we find that the autocorrela- 
tion times are roughly of the same order of magnitude, and may even become 
larger by multicanonical sampling. In fact, for O = m2 multicanonical sam- 
pling only increases the statistics in the exponentially suppressed tail of the 
canonical probability distribution P(m2). This observation, however, does 
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not affect our overall claim that the combination of multigrid techniques with 
multicanonical updating does significantly enhance the performance of the 
Monte Carlo process even though this improvement is practically irrelevant 
in the case of even observables. 

4.4 Real-time performance 

To conclude the analysis of the performance of the multicanonical multi- 
grid algorithm we finally need to look at the real-time work needed for 
the different algorithms. From a theoretical work estimate []TT| it follows 
that for 7 < 2*^ the additional work necessary to perform one complete W- 
cycle in comparison to a simple multicanonical sweep is given by a constant 
factor. For 7 = = 2 this factor is predicted to be close to 2. With 
our implementation on a CRAY X-MP we have measured updates times 
per site and cycle of Ueai = 11.2, 10.3, 9.5, Q.lytis for the W-cycle (W) resp. 
^reai = 4.0, 3.9, 3.9, 3.8/xs for the Metropolis (M) algorithm for lattices of size 
L = 8, 16, 32, 64. On a 128^-lattice our program would run with 8.2 fis (W) 
resp. 3.7 /is (M), and on a 256^-lattice with 8.1 /is (W) resp. 3.7 fis (M). 
It goes without saying that these numbers strongly depend on hardware fea- 
tures of the computer and on details of the implementation. We conclude 
that the gain in reduction of the autocorrelation times of a factor of 20 is 
roughly halved by the additional work needed to perform a W-cycle. Thus it 
is established that the combination with multigrid techniques enhances the 
performance of the standard multicanonical algorithm by about one order of 
magnitude, asymptotically independent of the linear lattice size L. 



5 Interface tension 

Having tested the performance of the algorithms we now turn to the evalua- 
tion of some observables of interest. Before doing so we recall that standard 
reweighting techniques |^ allow to compute expectation values of observ- 



ables for an appreciable range away from the simulation point. Since in 
the multicanonical case several different reweighting factors are employed we 
briefly review the histogramming technique for this case. 

In order to reweight to a new set of parameters we use the notation of 
eqs.(||,||) and notice that expectation values of canonical observables O = 
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0{Kq, Ml, M2, Mi) are obtained from the multicanonical distribution by 
computing 



;u;eani/x , Q) ^ dK^dM^dM^dM^ e/(^i/^)7V exp{-7i-,--} ^ ^ 

J dKodMidM2dMi P^'^^^Oe-f^^^/^^ 
J dKodMidM2dM4 p~ef(^^/^') ^ 

Here N — N{Ko, Mi, M2, M3) denotes the density of states for the variables 
Ko andM,, ^-.^(i^o, Mi, M2, M4) = Ko/2- {fj,y2)M2 + gM^ + f{Mi/V) is 
the multicanonical energy, and P^^g'^ oc A^exp{— H™2"^'^} is the multicanoni- 
cal probability distribution. Also we have dropped cancelling normalization 
factors. In a multicanonical Monte Carlo simulation configurations are sam- 
pled with a probability cx: Pj^^g^. Hence if we record the evolution series 
Mi of a simulation performed for one set of parameters {iJ>^,g) the expecta- 
tion value of an observable O for some other set of parameters {fi''^,g') can 
now in principle be calculated by multiplying with a reweighting factor. For 
example, the expectation value for /i'^ ^ fi^ would simply be given by 

2 ^ IdKodMidM2dM,P;^-^-Oef(^^/^)e'^^^ 

(0)can()U = — -72—2 ■ (32) 

/ dKodMidM2dM4 p~ef(^^^/^^e'^^^^ 

The only restriction for the reweighting procedure is given by the fact that 
the statistical accuracy of the data deteriorates if one reweights the data 
to a set of parameters far away from the simulation point. The problem is 
illustrated in Figs. 7(a-f). Fig. 7a shows the joint probability distribution 
P(mi, 771.2) for the multicanonical simulation at /i^ — 1.40 and L — 64, and 
Fig. 7b shows the same distribution after reweighting to the canonical case. 
These figures are directly comparable to Fig. 3c. Again we see in Fig. 7a the 
flat region between the peaks which allows the system to travel from states of 
negative to positive magnetization. Note that the histogram depicted in Fig. 
7a does produce the flat one-parameter histogram of Fig. 3c after integration 
over 1712- For the histogram reweighting. however, it is important to realize 
that also for the multicanonical situation of Fig. 7a a reweighting in the 
parameter /^^ shifts the histogram towards regions of smaller 1712 where the 
multicanonical statistics is as bad as a canonical simulation would be. After 
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reweighting to /x^ = 1.375 and 1.35 the resulting distributions are depicted 
in Fig. 7c resp. 7e for the multicanonical and in Fig. 7d resp. 7f for the 
canonical case. Comparing the multicanonical distributions of Figs. 7c and 
7e with the original smooth distribution of Fig. 7a one clearly sees that the 
histograms get increasingly noisy since the reweighting procedure suppresses 
the high statistics regions in Fig. 7a in favor of regions where only few 
configurations were sampled. Note that the normalization was adjusted in 
such a way as to show the peaks at same height. Consequently the ^-scale 
varies over many orders of magnitude in Figs. 7 (a-f), i.e., the maxima of the 
distributions vary as 4358 (7a), 3858 (7b), 11941 (7c), 282 (7d), 678 568 (7e), 
and 130 (7f). Although from Fig. 7e one would expect that the reweighting 
already breaks down it was nevertheless possible to find overlapping regions 
when reweighting our data in the intervals between the simulation points. 

Since we were mainly interested in the first-order phase transition we did 
not focus on thermodynamic properties at criticality. We only mention that 
the reweighting technique in principle allows us to compute the susceptibility 
X and specific heat C at the second-order transition line starting from our 
simulation data for = 1.30. In this way we have determined the transition 
point by extrapolating the finite lattice peak locations of x and C for L — oo 
and found a critical value of /i^ = 1.270(7) which is slightly larger than 
but still compatible with the value of /i^ = 1.265(5) found by Toral and 
Chakrabarti . 



More interesting in our context of an investigation of first-order phase 
transitions is a study of the interface tension. Here again we may use 
reweighting techniques. As discussed in section 2 the interface tension 
can easily be extracted from a histogram of mi by the relation 

1 pmax 

(^L = —r^^-^- (33) 

^ pram \ 

Since in the canonical distribution p™^^ is larger than p™™ by many orders 
of magnitude a reliable numerical evaluation of this relation is only possible 
for multicanonical simulations. In a canonical simulation there would only 
be very few configurations (if any) around P™™ and the relative statistical 
error of P™™ would be prohibitively large. Due to the flat multicanoni- 
cal histograms on the other hand the region around P™™ is sampled with 
the same statistical accuracy as the region around the maxima. A simple 
determination of the maximum resp. minimum of the histogram strongly 
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depends on the bin size of the histogram and tends to overestimate the ra- 
tio P'^^^l pmin for moderate bin sizes. To avoid this problem we determined 
pmax pmm |^-y g^^^jj^g parabolas to the extremal points of the histogram. 
For the histograms we used a bin size of 0.004, i.e., we had of the order of 
10^ entries in the bins between the maxima. For the fits of the maxima we 
cut the data at 0.85 x p^^^^ and for the fits of the minima we used data 
from mi = —0.2 . . . 0.2. We have checked that the results did not sensibly 
depend on the specific choice of the histogramming bin size or the cutting 
parameters for the fits. 

Fig. 8 shows the interface tension ctl for various lattice sizes L. The 
solid circles show the points where the actual simulations were performed, 
the interpolating lines were obtained by reweighting. Note that we have 
reweighted the data up to the mid points where the reweighted data from 
above meet those which were reweighted from below. Judging from Figs. 
7(a-d) we believe that this range still gives reliable values. For our small 
lattices L = 8, 16, and 32 we were also able to reweight our data well beyond 
the critical value for the infinite system. To obtain values for the infinite 
volume interface tension (Too we have extrapolated the (reweighted) data ac- 
cording to a fit of the form ctl = (Too + a/L. The squares in Fig. 8 show 
our infinite volume interface tensions at the simulation points. The precise 
values are listed in Table ^ 

From the universality with the two-dimensional Ising model we expect 
that the interface tension varies linearly with /x^ since for the Ising model the 
critical exponent u is equal to 1. Looking at the dependence of a^o with /i^ 
we find indeed that the interface tension a^o behaves like a^o = ax (fi"^ — fi^). 
A linear fit of the three data points of (Too intersects the /x^-axis at a value 

= 1.274(3) which agrees with our value obtained from extrapolating the 
maxima of the susceptibility x a-nd the specific heat C. For the interpretation 
of these data we would like to point out, however, that the goodness of the 
fit (Tl = (Too + a/L, which is perfectly satisfactory for large /x^, somewhat 
deteriorates as one approaches the critical line. In fact, for fi"^ = 1.30 a fit of 
the form (Tx, = aoo+a/ L+b/L^ gives a better chi-squared. Applying this fit to 
values of (Tl reweighted to values of /i^ larger than 1.31 on the other hand does 
not give a more consistent fit. The reason for this is probably the fact that 
for /i^ = 1.30 our histograms do not show a really flat region around mi ^ 
yet (cp. Fig. 2a). Hence interactions between the interfaces apparently are 
not yet completely negligible. It should also be kept in mind that in the 
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determination of (Too in the vicinity of /i^ quite a bit of numerical analysis is 
involved. Our extrapolation of the infinite volume interface tension doo to /i^ 
is therefore to be taken with a cautious mind. In particular, our data do not 
allow to decide how far away from the critical value /i^ the assumed linearity 
of (Too = a X (yU^ — ^1) actually holds. 



6 Concluding remarks 

We have shown that a combination of the multicanonical reweighting al- 
gorithm with multigrid update techniques reduces autocorrelation times of 
the Monte Carlo process at the field-driven first-order phase transitions of 
the two-dimensional 0^-model by a factor of ~ 20 when compared with stan- 
dard multicanonical Metropolis updating. Taking into account the additional 
work required for the multigrid W-cycle this effectively improves the real-time 
performance of the Monte Carlo process by about one order of magnitude 
compared with standard multicanonical simulations. 

Having established this gain in performance it would now be interesting 
to perform simulations of the (/)^-model in three or four dimensions as the 
immediate next step. Due to the generality of both the multicanonical for- 
mulation as well as the multigrid technique the algorithm is not restricted to 
only this one model and it is hoped that the method may further enhance 
Monte Carlo studies of first-order phase transitions or tunneling phenomena 
in quantum statistics [IG, F?]. 
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Appendix A: Error propagation for reweight- 

ing simulation data 

For any observable O (e.g. rrii = ^ Ya=i 4'i) expectation values in the canon- 
ical ensemble, (C)can5 are calculated as 

{0)c... = (Al) 

{w) 

where (. . .) (without subscript) denote expectation values with respect to the 
multicanonical distribution and w = exp(/) is the inverse reweighting fac- 
tor. ^ In a Monte Carlo simulation with a total number of measurements 
these values are, as usual, estimated by the mean values 

{Ow)^Ow = ^E^^^M (A2) 

Nm 

{w)^w = irpE^i' (A3) 

where Oi and Wi denote the measurements for the i-th configuration. Hence 
(C')can is estimated by 

(O)ean ^O^^. (A4) 
W 

The estimator O is biased, 



(d) = (OWl-f2^ + |;^ + ...l, (A5) 

{Ow){w) \w}{w) 

and fluctuates around (O) with variance, i.e. squared statistical error 



2 /.o\2 , {Ow;Ow) {w;w) ^ {Ow;w) 

€0 = ^ can + /-\2 +■■■ . A6 

^ {Ow)^ (w)^ {Ow){w) 



Here {Ow]w) = (Oww) — {Ow){w), etc. denote (connected) correlations of 
the mean values, which can be computed as 



9 int 



{Ow;w) = {0,wr,w,)^^, (A7) 



^Of course, the same considerations apply to the standard reweighting method p4[ as 
welL 
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where 

N, 



int _ int _ 1 V {OoWo;Wk) k 

-o.;. - r^;0. -2 + 1. ^Oowo; wo) nJ ^^^^ 

is the associated integrated autocorrelation time of measurements in the mul- 
ticanonical distribution. 

Hence the statistical error is given by 

2 ^ /.^v2 r {OiWu O^Wi) 2t'^^.(^^ {wf, Wj) 2t'^^ 

^ {OiWi]Wi) 2rg^.^ 

Since for uncorrelated measurements t^^.^^ = r^^.^ = r^°^ = 1/2 it is 
useful to define an effective multicanonical variance^ 



such that the error ( |A9|) can be written in the usual form 

4 - -f-^^ (All) 



with tq collecting the various autocorrelation times in an averaged sense. 
For a comparison with canonical simulations we need one further step since 

9 can 

= {^k)--^ (A12) 

but 0"^ 7^ (coi)'^^'^ — (C^iSC'j)- Hence we define an effective autocorrelation 
time through 

9_cff _cfl 

4 = (^J-'^f^ = i^'o)-^.^ (A13) 



I.e., 



r'^ = 77?L-o. (A14) 



^^2^^ Kan 



^ In the multicanonical distribution this is nothing but an abbreviation of the expression 
on the r.h.s. but not the variance in the multicanonical distribution. 
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For symmetric distributions and odd observables we have {OiWi) = and 
this simphfies to 



such that 



ro = T^l,o^, (A16) 

and 

"^o = 7;;t^Jow;ow (A17) 
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Table Captions 



Tab. 1: Canonical simulation: Integrated and exponential autocorrelation 
times r™* and r*^^P, and flipping time t^'p using the standard Metropolis 
algorithm (M) or the multigrid W-cycle (W), /x^ = 1.30. 

Tab. 2: Multicanonical simulation: Integrated and exponential autocorre- 
lation times r'"* and t^^^ using the standard Metropolis algorithm (M) 
or the multigrid W-cycle (W). 

Tab. 3: Multicanonical simulation: Integrated and exponential autocorre- 
lation times T^l^f and T^^^f using the standard Metropolis (M) or the 
multigrid W-cycle (W). Also listed are the effective multicanonical vari- 
ance 0"^ and the canonical variance {cq.Y^^ for O = mi. From these 
the effective autocorrelation time for the canonical statistical error 
estimate can be computed according to eq.(^). For comparison, in 
the last column we also list the "flipping" time T^f for the diffusion 
between the peak maxima. Same values of ng as in Table 

Tab. 4: Multicanonical simulation: Mean values, variances, and covariance, 
as well as integrated autocorrelation times t^^.^^, r™^, and t^^.^ for 
O = 1712 and w = exp(/) and /i^ = 1.30. These values enter the statis- 
tical error estimate eq.(p5D for the even observable m2 which allows to 
compute the squared canonical error estimates and the corresponding 
effective autocorrelation time deflned in eq. (p^ . For comparison 
the same quantities were also obtained by direct jackkniving (jack). 

Tab. 5: Interface tension for various lattice sizes and /i^ = 1.30, 1.35, 
and 1.40. The inflnite volume interface tension (Too was obtained by a 
flt according to = a^o + a/L. 
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= mi 


= m2 


L 






int 
' mi 


exp 
' mi 


' mi 


He 


int 
'm2 


exp 
' m2 


4 


M 


3 


154.6(3.7) 


159(11) 


200.5(2.4) 


3 


12.366(82) 


15.07(34) 


4 


W 


1 


15.15(19) 


15.43(36) 


14.49(14) 


1 


1.4794(87) 


2.549(82) 


8 


M 


20 


847(19) 


877(54) 


1007(11) 


20 


28.19(15) 


40.2(1.3) 


8 


W 


1 


40.20(81) 


39.8(21) 


48.24(48) 


1 


1.841(16) 


3.70(14) 


16 


M 


150 


5780(110) 


5710(320) 


6575(62) 


2 


55.4(2.3) 


118(14) 


16 


W 


8 


239.6(3.8) 


248(12) 


275.9(2.3) 


1 


3.475(43) 


8.06(27) 



Table 1: Canonical simulation: Integrated and exponential autocorrelation 
times r™* and t^'^p, and flipping time r^f using the standard Metropolis 
algorithm (M) or the multigrid W-cycle (W), //^ = 1.30. 
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3400(270) 


3000(630) 


194(75) 


820(780) 



Table 2: Multicanonical simulation: Integrated and exponential autocorrela- 
tion times r"^* and r^^^ using the standard Metropolis algorithm (M) or the 
multigrid W-cycle ( W) . 
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6.90(22) 


0.6275(82) 
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Table 3: Multicanonical simulation: Integrated and exponential autocorrela- 
tion times T™*g/ and T^^^f using the standard Metropolis (M) or the multigrid 
W-cycle (W). Also listed are the effective multicanonical variance and the 
canonical variance {o^oj'^^'^ for C> = mi. From these the effective autocorre- 
lation time for the canonical statistical error estimate can be computed 
according to eg. (^61) . For comparison, in the last column we also list the 
"flipping" time r^f for the diffusion between the peak maxima. Same values 
of Tie as in Table 0. 
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Table 4: Multicanonical simulation: Mean values, variances, and covari- 
ance, as well as integrated autocorrelation times t^^.^^, t™^, and t^^.^ for 
O = 777.2 and w = exp(/) and /i^ = 1.30. These values enter the statistical 
error estimate eq.(p5D for the even observable 7772 which allows to compute 
the squared canonical error estimates and the corresponding effective au- 
tocorrelation time defined in eq . (p4|) . For comparison the same quantities 
were also obtained by direct jackkniving (jack). 
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L 


fi' = 1.30 


= 1.35 


^JL^ = 1.40 


8 


0.14826(58) 


0.19013(52) 


0.23668(79) 


16 


0.10526(47) 


0.15634(51) 


0.21288(64) 


32 


0.07095(39) 


0.12690(42) 


0.18964(49) 


64 


0.05173(37) 


0.11260(50) 


0.17732(61) 


oo 


0.03443(47) 


0.09785(60) 


0.16577(73) 



Table 5: Interface tension ul for various lattice sizes and /x^ = 1.30, 1.35, 
and 1.40. The infinite volume interface tension a^o was obtained by a fit 
according to ctl = + a/ L. 
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Figure Captions 



Fig. 1: Interfaces in a typical mixed phase configuration for /x^ = 1.40 and 
L = 64. Values of (pi > 0.5(< —0.5) have been depicted black (white), 
values of < 0.5 have been depicted by varying grey shades. For 
this configuration the magnetization was mi ^ —0.2. 

Fig. 2a-c: Flat multicanonical distributions as compared to the canonical 
double-peak histograms of mi after reweighting for //^ = 1.30 and dif- 
ferent lattice sizes L = 8, 16, and 32. 

Fig. 3a-c: Flat multicanonical distributions as compared to the canonical 
double-peak histograms of mi after reweighting for L = 64 and /i^ = 
1.30, 1.35, and 1.40. The arrow in Fig. 3c indicates the value of mi 
which was measured for the configuration displayed in Fig. 1. 

Fig. 4a-d: Evolution series of mi for L = 16 and /x^ = 1.30 using the 
canonical Metropolis (Fig. 4a), the multicanonical Metropolis (Fig. 
4b), the canonical multigrid W-cycle (Fig. 4c), and the multicanonical 
multigrid W-cycle (Fig. 4d). In Figs. 4a and 4c we also show the 
digitalized time series according to the simple two-state flip model. 
The time scales were adjusted so that an evolution over roughly 30r™* 
is displayed in each figure. 

Fig. 5: Autocorrelation time r(A;) together with a three-parameter fit ac- 
cording to eq. (|T^) for O = miexp(/), /i^ = 1.40, L = 64, and the 
multicanonical W-cycle. The horizontal line shows the integrated au- 
tocorrelation time T™*g/ = r(oo) together with error bounds. The inset 
shows a logarithmic plot of the autocorrelation function A{j) together 
with linear fits \nA{j) = const — j/'T^^^.f- See the text for a detailed 
explanation. 

Fig. 6: Effective autocorrelation times and integrated autocorrelation 
times T™*g/ as a function of L. Straight hues are fits according to 
r = aL^. 
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Fig. 7a-f: Two-dimensional histograms of mi and 1112 simulated for L = 64 
and /i^ = 1.40 and reweighted to different values of /x^. The multi- 
canonical distributions are shown in Fig. 7a for //^ = 1.40 without 
reweighting, in Fig. 7c after reweighting to /i^ = 1.375, and in Fig. 
7e after reweighting to fi^ = 1.35. Figs. 7b, 7d, and 7f show the re- 
spective canonical distributions after additional reweighting with the 
multicanonical reweighting factor exp(/(m)). 

Fig. 8: Interface tension (Tl as a function of fi^ for L = 8, 16, 32, and 64. 
The filled circles show the actual simulation data and the dashed lines 
were obtained by reweighting. The values for a = cJoo were obtained by 
an extrapolation according to aL — (7 00 -\- a/L, and the solid straight 
hne shows a fit (Too = a{iJ? — nl) with nl — 1.274(3). 
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Figure 1: Interfaces in a typical mixed phase configuration for fi"^ = 1.40 
and L = 64. Values of 0j > 0.5(< —0.5) have been depicted black (white), 
values of |0j| < 0.5 have been depicted by varying grey shades. For this 
configuration the magnetization was mi ~ —0.2. 
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Figure 2: Flat multicanonical distributions as compared to the canonical 
double-peak histograms of mi after reweighting for //^ = 1.30 and different 
lattice sizes L = 8, 16, and 32. 
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Figure 3: Flat multicanonical distributions as compared to the canoni- 
cal double-peak histograms of mi after reweighting for L = 64 and = 
1.30, 1.35, and 1.40. The arrow in Fig. 3c indicates the value of mi which 
was measured for the configuration displayed in Fig. 1. 

45 




50000 100000 150000 




5000 10000 15000 




2000 4000 6000 8000 




200 400 600 800 1000 



Figure 4: Evolution series of mi for L — 16 and //^ — 1.30 using the canonical 
Metropolis (Fig. 4a), the multicanonical Metropolis (Fig. 4b), the canonical 
multigrid W-cycle (Fig. 4c), and the multicanonical multigrid W-cyclc (Fig. 
4d). In Figs. 4a and 4c wc also show the digitalizcd time series according 
to the simple two-state flip model. The time scales were adjusted so that an 
evolution over roughly 30t™* is displayed in each figure. 
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Figure 5: Autocorrelation time T{k) together with a three-parameter fit 
according to eq. (16) for O = mi exp(/), = 1.40, L = 64, and the multi- 
canonical W-cycle. The horizontal line shows the integrated autocorrelation 
time T^l^f = t{oo) together with error bounds. The inset shows a loga- 
rithmic plot of the autocorrelation function A{j) together with linear fits 
In A{j) — const — j/T^^^f- See the text for a detailed explanation. 
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Figure 6: Effective autocorrelation times and integrated autocorrelation 
times T^jg/ as a function of L. Straight lines are fits according to r = . 
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Figure 7: Two-dimensional histograms of mi and m2 simulated for L = 64 
and /i^ = 1.40 and reweighted to different values of fi"^. The multicanonical 
distributions arc shown in Fig. 7a for /x^ = 1.40 without reweighting, in Fig. 
7c after reweighting to //^ = 1.375, and in Fig. 7e after reweighting to /x^ = 
1.35. Figs. 7b, 7d, and 7f show the respective canonical distributions after 
additional reweighting with the multicanonical reweighting factor exp(/(m)). 
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Figure 8: Interface tension gl as a function of //^ for L = 8, 16, 32, and 
64. The filled circles show the actual simulation data and the dashed lines 
were obtained by reweighting. The values for a = (Too were obtained by an 
extrapolation according to (Tl = + a/L, and the solid straight line shows 
a fit (Too = a(/x2 - fxl) with /x^ ^ 1.274(3). 
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